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ABSTRACT 

We present sensitive, arcsecond-resolution Submillimeter Array observations of the 
O J ^^CO J=2— 1 line emission from the circumstehar disk orbiting the double-hned spec- 

troscopic binary star V4046 Sgr. Based on a simple model of the disk structure, we 
use a novel Monte Carlo Markov Chain technique to extract the Keplerian velocity 
field of the disk from these data and estimate the total mass of the central binary. 
Assuming the distance inferred from kinematic parallax measurements in the literature 
(d ^ 73 pc), we determine a total stellar mass = 1-751q Qg M© and a disk inclination 
id = 33?5^5'4 from face-on. These measurements are in excellent agreement with inde- 
■ pendent dynamical constraints made from multi-epoch monitoring of the stellar radial 

velocities, confirming the absolute accuracy of this precise (^few percent uncertainties) 
disk-based method for estimating stellar masses and reaffirming previous assertions that 
the disk and binary orbital planes are well aligned (with | 0.1±1°). Using these 

results as a reference, we demonstrate that various pre-main sequence evolution models 
make consistent and accurate predictions for the masses of the individual components of 
the binary, and uniformly imply an advanced age of ^5-30 Myr. Taken together, these 
results verify that V4046 Sgr is one of the precious few nearby and relatively evolved 
pre-main sequence systems that still hosts a gas-rich accretion disk. 

Subject headings: protoplanetary disks — stars: individual (V4046 Sgr) 



1. Introduction 

Mass is the fundamental property that sets the evolutionary path of a star. The masses of 
young stars are of particular interest in many astrophysical problems: they provide unique in- 
formation about th e star formation process (e.g., accretion histories, the initial mass function; see 



Bastian et al.ll2010l ) and are thought to have a substantial influence on the evolution of their circum 



stella r material, and therefore the efficiency of the planet formation process (e.g., see lAlibert et al. 



201 ll ). Unfortunately, measurements of pre-main sequence (pre-MS) star masses are difficult and 



accordingly rare. Unlike their more evolved counterparts, the location o f a young star in the 



Hertzsprung- Russell (HR) diagram does not provide a robust estimate of (IHillenbrand fc White 
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2004 ) . The theoretical mod els for stehar evolution at these early st ages are pla gued with uncer- 



tainties related to rotation (ISiess fc Livid 119971 : 



BarafFe et al 



2009 



BarafFe fc Chabrier 



Mendes et al 



19991). accretion (jSiess et alJ 119971 : 



2010l). magnetic fie lds (jD'Antona et alJl2000l ). atmosphere 



properties (e.g., convection, opacities; iBarafFe et al.ll2002l ). and unknown initial conditions. Ul- 
timately, a more nuanced understanding oF star and planet Formation requires pre-MS evolution 
models that are empirically calibrated with direct, independent, and accurate measurements. 

The only direct methods available For measuring stellar masses are based on orbital dynamics. 
In sufficiently close pre-MS binary star systems, can be estimated From the stehar o rbits us- 
ing multi-epoch radial velocity (R V) measurements (e.g .. 



long-term astrometric monitoring ( iTamazian et al.ll2002l : ISchaeFer et al. 



Mathieu et a" 



20031, M 



19971) and /or 



Duchene et al. 



20061 ) . For a double-hned spectroscopic binary (SB2), the RV method provides a robust estimate 
oF M*(sinz*)^ For each component, where is the inclination angle oF the orbit projected on the 
sky. For "visual" binaries, the astrometric monitoring technique ofFers a constraint on the quantity 
Mtotd^, where Mtot is the sum oF the stellar masses and d is the distance to the binary. Gener- 
ally, the stellar mass estimates From these methods are inherently uncertain due to their strong 
dependences on the unknown values oF or d. The {M*, i^} degeneracy is broken For the spe- 
cial c ase oF an eclipsing; SB2, but Few pre-MS systems with such Favorable orientations are known 
(e.g.. lStassun et al.ll2004l : iMorales-Calderon et al.ll2012L and reFerences therein). In a subset oF ideal 
cases, the RV and astrometric techniques can be combined to aheviate the uncertainties related to 



U*, d} and extract acc urate values (e.g., IStefFen et al.ll200ll : ISchaeFer et al.ll2008l : iBoden et al. 
2OO5I . I2OO7I . I2OO9I . l2012h . 



These standard methods are only applicable For binary stars with a narrow range oF orbital 
separations. Alternatively, For any isolated young star with a circumstellar disk, M^^ can be de- 
termined From a single millimeter-wave interFerometric observation oF an optically thick emission 
line (with a linear dependence on d). This latter technique relies on modeling the spatially and 
spectrally resolved Keplerian rotation curve oF the molecular ga s disk that orbits the young star 
(IKoerner et al.lll993l : IPutrev et al.lll994l . Il998l : ISimon et al.ll2000l ). While this method has extraor- 
dinary value in its more general applicability, it has only been successFully employed For small 
samples. Attempts to expand its reach have been Frustrated by molecular cloud contamination 
and observational limitations in resolution and sensitivity. Moreover, a reconstruction oF the disk 
velocity fi eld necessarily involves fitt ing; such data with a relatively complicated model oF the disk 
structure (IBeckwith fc Sargentlll993l ). Given that added complexity, there i s naturally sonie con - 



cern about the absolute accuracy oF the estimates From this method (e.g. . iGennaro et al. l2012h 



despite the impressive Formal precision oF the measurements (~2-3%; e.g.. IPietu et al. 2007). 



The young binary V4046 Sgr provides a rare opportunity to benchmark the disk kinematics 
method For estimating against the more traditional RV technique For a SB2 system. V4046 Sgr is 
a nearly equal mass (q ^ 0.94) pair oF solar- type pre-MS stars in a circular (e < 0.01), non-eclipsing 



orbit with a 2.4 day period (a sm ^ 5.1 R0; iBvrnd 119861 : iQuast et al.ll2000l : IStempels fc Gahml 



20041 ) . The system is completely isolated From any known molecular clouds and has been kinemati- 
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cally associated with the ^8- 20 Myr-old Pic m oving; group; a moving cluster analysis suggests it is 
relatively nearby, d = 73pc ( Torres et al. 20061 ) . Despite its advanced age, V4046 Sg r hosts a large 



and massive circumbinary disk that exhibits a rich molecular emission line spectrum (iKastner et al 



20081 : [Rodriguez et alJl2010l : lOberg et al.ll201ll ) . Since the binary orbit is tigh t and circular, it has no 



dyna mical impact on the disk structure outside a radius of ^0.1 AU (e.g., see I Artymowicz fc Lubow 



1994 ). With its central SB2 host, rare proximity to the Sun, lack of molecular cloud contamination, 
and intrinsically bright, spatially extended line emission, the V4046 Sgr disk is an ideal target to 
assess the accuracy of the disk kinematics technique for measuring M*. 



In this article, we build on some initial work by [Rodriguez et al.l (120101 ) and use high-quality 
spatially and spectrally resolved observations of the CO J=2— 1 emission line to measure the 
velocity field of the V4046 Sgr disk and extract the total mass of the close binary at its center. Our 
millimeter-wave observations with the Submillimeter Array (SMA) and data calibration procedures 
are described in 321 A detailed overview of the modeling analysis is provided in 331 The modeling 
results ar e pres ented and compared with the complementary RV analysis of the central SB2 by 
Stempelsl (120121 ) in 311 These results are discussed in the context of the V4046 Sgr system in 
particular, pre-MS evolution models more generally, and future prospects for estimates from 
the disk kinematics method in 33 Some key conclusions from this work are summarized in 321 



2. Observations and Data Reduction 

The V40 46 S^i circumb inary disk was observed at 225 GHz (1.3 mm) with the Submillimeter 



Array (SMA: lHo et al.ll2004l ) on four occasions starting in 2009, using each of the available antenna 
configurations: sub-compact (baseline lengths of 9-25 m; 2011 Mar 18), compact (16-70 m; 2009 
Apr 25), extended (28-226 m; 2009 Feb 23), and very extended (68-509 m; 2011 Sep 4). The 
SMA double sideband receivers were tuned to simultaneously observe the J=2— 1 transitions of 
i^CO, i^CO, and C^^O at 230.538, 220.399, and 219.560 GHz, respectively, and the adjacent dust 
continuum. The correlator was configured to place those emission lines in separate 104 MHz spectral 
chunks and sample them finely with 512 channels per chunk, corresponding to a native velocity 
resolution of ^0.25 km s~^. The continuum was observed with a more coarse frequency sampling (in 
3.25 MHz channels), with a total bandwidth of 1.6 and 3.6 GHz in 2009 and 2011, respectively. The 
observations cycled between V4046 Sgr and the quasars J1924-292 (15° away) and J1733-130 (22° 
from V4046 Sgr, 30° from J1924-292) with a total loop time of 10-20 minutes. The bright quasars 
3C 84 and 3C 454.3 were also observed as bandpass cahbrators, along with Uranus, Callisto, and 
Ceres for use in determining the absolute scaling of the amplitudes. All of the data were collected 
in outstanding weather conditions for this observing frequency, with an atmospheric zenith optical 
depth of only ^0.05 (corresponding to a precipitable water vapor level of ^1.0 mm). 

Each individual dataset was calibrated independently with the MIR software package. The 
bandpass response was corrected based on observations of bright quasars, and broadband continuum 
channels were generated from the central portions of all line-free spectral chunks. The visibility 
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Fig. 1. — Naturally weighted channel maps of CO J=2— 1 emission from the V4046 Sgr disk. The 
synthesized beam dimensions are marked in the bottom left. In each 0.4 km s~^-wide channel (LSR 
velocities labeled in the lower right corner of each panel), contours are drawn at 0.21 Jy beam~^ 
(3(j) intervals, starting at 0.18 Jy beam~^. A linear color scale represents the mean velocity of each 
channel, with red and blue denoting receding and approaching velocities, respectively. 

amplitude scale was derived by bootstrapping the gain calibrator (quasar) flux densities from the 
observations of Uranus, Callisto, or Ceres, with a systematic uncertainty estimated at 10-15%. 
The antenna-based complex gain response of the system was determined with reference to J 1924- 
292, and the quality of the phase transfer was assessed using the observations of J 1733- 130. That 
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Fig. 2. — (left) The CO J=2— 1 moment maps for the V4046 Sgr disk, on the same angular scale 
as in Figured! The 0^^ moment (velocity-integrated intensity) map is overlaid in contours shown 
at 5(7 intervals (0.35 Jy km s"-*^ beam"-*^), starting at the 3 a level (0.21Jy km s"-*^ beam"-*^). The 
-yst jxioment (intensity- weighted velocities) map is shown in color, with a scale bar for reference. 
{right) The integrated CO J=2— 1 line profile from the V4046 Sgr disk, constructed from a square 
box 10'^ on a side centered at the stellar position, and assuming a 2 a threshold for inclusion. 



comparison suggests only a small amount of "seeing" (^O'.'l) was introduced by atmospheric phase 
noise (or small baseline errors), consistent with the excellent observing conditions. After applying 
the appropriate (and small) p hase shifts to accoun t for the V4046 Sgr proper motion (/Xq, cos S = 
0^.^003 yr~^, ij^s = —0^.^052 vr~^: IZacharias et al.ll2010l ) and confirming that the continuum amplitudes 
from different array configurations were consistent on overlapping baseline lengths, the visibility 
datasets from each observation were combined. The observations of the dust continuum and CO 
isotopologue emission will be presented in a separate article; the focus here will be s olely on the ^^CO 
J=2— 1 emission. Note that although the 2009 data were originally presented by [Rodriguez et al. 



fcoiol ). those data have been re-calibrated here for consistency (and modest improvements). 

The CO visibilities were continuum-subtracted and truncated outside a projected baseline 
length of 200 kA to reduce the data volume and improve the signal-to-noise ratio. They were then 
Fourier inverted, deconvolved with the CLEAN algorithm, and restored with a synthesized beam 
using the MIRIAD package. The naturally-weighted spectral images shown as channel maps in 
Figure 1 were synthesized on a 0.4 km s~^ smoothed velocity scale with a 1^.^55 x 1^.^29 beam (at 
a position angle of 24°). The typical RMS noise level in each channel is 70mJy beam~^. There 
is CO emission firmly detected (>3a) out to ±4.8 km s~^ from the systemic velocity, estimated 
to be i;(LSR) = +2.87 =b 0.05 km s"-*^ (corresponding to —6.26 =b 0.05 km in the heliocentric 
frame), with an integrated intensity of 27.7=b2.8 Jy km s~^ and a peak flux density of 1.48=b0.16 Jy 
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beam~^ (17±2K; a peak S/N = 21), including the calibration uncertainties. Figure [2] shows a 
map of the velocity-integrated CO intensities (0*^ moment; contours) overlaid on the intensity- 
weighted velocities (1^* moment; colorscale)^ as well as a spatially integrated CO spectrum. These 
data exhibit a molecular gas disk with a clear rotation pattern, from east (blueshifted) to west 
(redshifted) , and suggest a modest inclination angle to t he line of sight. Near t he systemic velocity, 
the CO emission subtends ^5^' in radius (^365 AU; see iRodriguez et all 12010 ). 



3. Modeling Analysis 



The fundamental goal here is to derive a dynamical estimate of the central stellar mass based 
on the kinematic properties of the V4046 Sgr gas disk. In order to extract from a measurement 
of the disk velocity field, we need to constr uct a detailed physical mod el of the disk structured 
We adopt a modeling formalism motivated by lBeckwith fc SargentI (|l993l ), which makes three basic 
assumptions about disk properties. First, the disk material is assumed to be orbiting in Keplerian 
rotation around a point mass, meaning the central stars are treated as one object that dominates 
the disk velocity field. Previous work has found little evidence to support (simple parametric) 
deviations from Keplerian i;-fields in disks (e.g., see Simon et al. 2000); the corollary criterion that 
the disk mass is only a small fraction of the stellar mass {Md/M^ <C 1) can be confirmed a posteriori 



(see S 4). Second, the disk is assumed to be geometrically thin at all radii. Although IPietu et al 



(120071 ) argued that this may not be valid at very large radii (>800AU), the higher for the 
V4046 Sgr binary (which decreases the disk scale height; see below) and the limited extent of the 
CO emission from its disk (<400AU) suggest it is a reasonable approximation in this case. And 
third, in the context of the disk structure, we assume that the gas is vertically isothermal and in 
hydr ostatic pressure equili brium. A more realistic model would include a temperature inversion 
(e.g.. IPAlessio et al.lll998l ). However, excluding that kind of added complexity does not diminish 
the accuracy of an M^, estimate: the emission from a single CO line is g:ene rated in a narrow ve rtical 
layer of the disk atmosphere, which has a roughly constant temperature (iDartois et al.ll2003l ). 



The model for the gas density distribution is constructed in cylindrical coordinates (r. 



p(r, z) 



exp 



2 \Hr 



(1) 



where S is the radial surface density profile and Hp is the pressure scale height at each radius. We 
assume a parametric versio n of the former that is appropriate for an accretion d isk with a static, 
power-law viscosity profile (jLynden-Bell fc Pringld Il974l : iHartmann et al.l 1 19981 ) and is currently 



^For the sake of convenient and general notation, we will refer to the central stellar mass as M*. In this particular 
case, M* corresponds to the sum of the stellar masses in the V4046 Sgr binary. 
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the basis for most dust-based disk density measurements (e.g.. I Andrews et alJl2009l . l2010l ). 



exp 



2-7' 



(2) 



where 7 is a density gradient, Tc is a characteristic radius, and Sc is a normahzation equivalent to 
e • S(rc). The scale heights are calculated with the explicit assumption of hydrostatic equilibrium 
and a constant vertical temperature profile, such that 



kT 



1/2 



(3) 



where is the sound speed, is the angular velocity, T is the radial temperature profile, k is 
the Boltzmann constant, G is the gravitational constant, mn is the mass of a hydrogen atom, and 
fi — 2.37 is the mean molecular weight of the gas. We further assume a simple power-law behavior 
for the radial temperature profile, 

where g is a temperature gradient and Tio is the gas temperature at a radius of 10 AU. The 
parametric descriptions in Equations (l)-(4) completely characterize the temperature and density 
structure of a model gas disk. The gas kinematics are described by a Keplerian velocity field. 



GM^ 



Vr v. 



0, 



(5) 



meaning there is only rotation, and no net motion in the radial or vertical dimensions. 



Given a physical disk structure, we can construct a model of the CO J=2— 1 emission with 
the added assumption that the energy levels relevant for this transition are populated according 
to local thermodynamic equilibrium (LTE). A l thoug h the LTE approximation is not always valid 
in protoplanetary disks, iPavlyuchenkov et al.l (120071 ) demonstrated that it is an appropriate sim- 
plification for the low energy and high optical depths associated with this particular species and 
transition. We model the emission line intensity distribution by integrating the radiative transfer 
equation along each sight line 5, so that 



roo 

Jo 



(6) 



where Ky{s) is the absorption coefficient, Sy{s) — By{T) is the source function (here equivalent 
to the Planck function), and Ty{s) = Kj^{x)dx is the optical depth along the line-of-sight. The 
absorption coefficient is the sum of contributions from dust and gas. For the former, we assume 
i^jy(5)dust — p{s) where the dust-to-gas mass ratio is = 0.01 and the grain opacity is /^^^ ^ 
(z//100GHz)cm2 g-i ( Beckwith et al.lll99o l). In the particular case of the V4046 Sgr disk, the 
contribution of dust to the absorption coefficient is effectively negligible compared to the gas; the 
specific choices of ( and hCj^ (within reasonable limitations) have no tangible effect on our results 
(moreover, continuum emission is removed from the data before our modeling analysis). 
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The absorption coefficient of the CO gas is calculated from the transition cross section weighted 
by the population of the lower energy level, ^, such that Kjy{s)co = n£{s)ajy{s) (note that in this 
specific case, ^ = 1). Since we assume the line is thermally populated in LTE, the level populations 
are determined by the local disk temperature via the Boltzmann equation. 



Ep 



(7) 



where is the transition energy, = 2^ + 1 is the statistical weight, Z is the partition function, 
and Xco is the CO fractional abundance (assumed to be constant everywhere in the disk). The 
details of the emission line model are encoded in the absorption cross section. 



where is the line profile function and the integrated cross section is 

^2 



(Jo = -— • n2i 

47r Qp 



9i+i 



87rz/2 



A 



21, 



(8) 



(9) 



where we have used the Einstein relation in the last equality to express the Einstein -^ coefficient 



in te rms of the Einstein-A coefficient provided by the LAMDA molecular database (ISchoier et al 



20051 ) . The line profile function naturally determines the shape of the emission line, where a given 
frequency jy corresponds to an intrinsic Doppler velocity = (c/z^o)(^ — ^o) relative to the line 
center, uq. The gas in a given disk model has a projected, line-of-sight velocity field given by 
'^obs ^ '^(/)(^) cos(/)sinZ(^, where is the inclination of the disk relative to the observer (such that 
id = 90° is edge-on; for the geometry, see llsella et alJl2007l ). The line profile shape at each frequency 
is determined by the difference between the Doppler and line-of-sight velocities. 



Ms) 



^/7TUoAv 



exp 



Av 



(10) 



and the effective line width, Av. The latter is comprised of the quadrature sum of contributions 
from thermal and non-thermal broadening terms. 



A., = 



rrir. 



1/2 



(11) 



where ruco is the mean molecular weight of CO and ^ is the contribution from microturbulence. 
We assume that the turbulent velocity width ^ is constant throughout the disk. 



The iBeckwith fc SargentI (119931 ) model formalism highlighted above is admittedly complex. 
A single model is completely described by 14 parameters: four quantify the distribution of CO 
densities {Xco, Sc, Rc^ 7}, two characterize the temperature structure (and therefore the vertical 
density structure) {Tio, q}, three play pivotal roles describing the disk kinematics {M*, ^, z/q}, and 
the remaining five relate the model to the observations - i.e., convert from physical coordinates in 
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the disk to the observed plane projected on the sky - including the disk inclination and major axis 
position angle {z^, PA^}, distance {d}, and disk center {ao, 60}. Some of these parameters can be 
reliably determined in a simple way, and then fixed to facilitate the data modeling problem (with 
no tangible impact on the quality of the determination). In this case, the systemic LSR velocity 
(effectively uq) was measured directly from the SMA data and set to vq = +2.87 km s"-*^. The disk 
center coordinates were estimated from both the dust continuum and CO emission morphology 
at the systemic velocity (see Fig. [1]) and set to ao = 18^14^1(^48, 60 = -32°4 r35:^08 (J2000), 



coincident with the composite V4046 Sgr stellar position in the UCAC3 catalog (jZacharias et al 



horn ). We adopted the distance inferred by lTorres et al.l (j2006l ) from the kinematic parallax (moving 



cluster) technique, d = 73pc (but see §5 regarding alternative values). 

Some additional practical simplications can be made, since the focus here is on the velocity field 
and not the disk structure details. Because the line opacity is so much greater than the dust opacity, 
the normalization parameters Xco and are not independent: we can only constrain their product. 
In practice, we adopt a joint CO disk mass parameter, Mco = XcoSc(27rr^) / (2 — 7), for that purpose. 
Moreover, after extensive experimentation, we concluded that two more parameters should be fixed. 
First, we found that the disk orientation was very well-determined (within 1°) at PA^^ = 76° (E 
of N), so that continued iteration on a more precise value was a waste of computational resources: 
fixing this parameter has negligible quantitative impact on the other model parameters or their 
uncertainties. And second, we fixed the surface density gradient, 7=1. The key parameters that set 
the density profile, {7, rd are anti-correlated and strongly degenerate; statistically indistinguishable 
model fits can be found over a large range of these parameters. The degeneracy is remarkably 
narrow, and could easily be missed with standard minimization algorithms (resulting in local 
minima with misleadingly tight constraints on the gradient). Fortunately, this degeneracy has 
minimal impact on the parameters most relevant for characterizing the disk velocity field: the 
precision and accuracy of a determination are not notably affected by the uncertainties in {7, 
Tc}. With that in mind, we have made the modeling process more tractable with a fixed 7. 

Making use of those simplifications, a synthetic CO spectral datacube can be calculated by 
specifying 7 free parameters, {Mco, ^c? ^lo? ^ci? M^}. That model datacube is then resampled 
at the observed velocities and spatial frequencies and processed in the same way as the data (see §2) 
to produce a set of synthetic spectral visibilities. Those model visibilities are evaluated with respect 
to the data by computing a composite value, summed over the real and imaginary components 
in 25 spectral channels. The best-fit set of model parameters and their associated uncertainties 
were determined with a Monte Carlo Markov Chain (MCMC) technique, using the affine invariant 



ensemble sampler developed by iGoodman fc Weard (I2OIOI . see also Foreman-Mackey et al. 2012), 



using a jump probability oc e~^^^/^. To our knowledge, this is the first application of these more 
sophisticated MCMC methods for parameter estimation in this particula r context: previous studies 
have relied on downhill simplex routines (e.g., Guilloteau fc Dutrey 19981: Isimon et ah 2000 ). some- 



times with clever modifications to address asymmetric uncertainties (Pietu et al...20(y^ )~ Although 
comparatively the MCMC technique is computationally expensive, this Bayesian treatment has the 
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Fig. 3. — Channel maps of the CO J=2— 1 emission hne from the SMA observations (top), the 
best-fit disk model (top), and the imaged residuals (bottom) with the channel velocity increasing 
from the top left to the bottom right. The last panel shows the 0*^ (contours) and 1^* (colorscale) 
moment maps as in Figure [2l The best-fit model parameters are listed in Table [H 

distinct advantage of providing the posterior probability distributions for each parameter in the 
complex, multi-dimensional parameter-space of the underlying disk model. 



4. Results 

The best-fit model derived from the SMA data is compiled in Table [H which lists the mode 
(peak) of the posterior probability distribution for each parameter along with its uncertainty, 
based on the range of values that encompass 68% of the distribution area (i.e., la uncertain- 
ties for a Gaussian distribution). Figure [3] makes a direct comparison between the data and the 
best-fit model in the spectral image plane, demonstrating the fit quality with only low-level (sta- 
tistically insignificant) residuals present in the channel maps. The best-fit model has a reduced 
= 967, 427/977, 000 = 0.9902. A more detailed view of the multi-dimensional results of the mod- 
eling analysis is presented in Figure [H where we have taken the 30,000 MCMC samples computed 
and marginalized over parameter subsets to display one and two-parameter posterior probability 
distributions. The marginalized distributions for individual parameters are single-peaked, while 
the paired distributions highlight internal degeneracies in the model. For example, there are tight 
correlations between normalizations and sizes (e.g., {Mco, rd, to conserve the integrated line in- 
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tensity; this would be modified were 7 a free parameter) or gradients (e.g., {Tm, gj, to main tain 



an appropriate temperature in the outer disk; iMundy et al.lll996l : I Andrews fc Wihiamsl 120071 ) . as 
weh a s more subtle associations amongst parameters related to line broadening (amongst rc, Tio, 
and ^: ICuilloteau fc Dutrevlll998l : iHughes et allboill ). The degeneracy of most interest here is the 
{M*, id} anti-correlation, which is a natural consequence of reproducing the observed line-of-sight 
velocity pattern, Vqi^s oc Vksinid oc \^M^sinid (i.e., oc 1/ sin^ id). 

The best-fit model parameters that describe the V4046 Sgr disk structure are typical for young 
protoplanetary disks. The adopted (fixed) surface density gradient (7 = 1) and characteristic 
radius ( tc. = 45^3 AU) lie ne a r the median values for a survey of disk structures in the Ophiuchus 
region JAndrews et all bood . I2OI0I ): the temperature profile (Tio = 115 ± 5K, g = 0.63 zb 0.01; 
and therefore vertical density distribution) and turbulent linewidth (^ = 0. 14 zb 0.01 km s~^) are 
comparable to those inferred from CO eniission in other T Tauri disks (e.g., iGuilloteau fc Dutrey 



199 



If we assume the CO mass fraction found in dark 



Pietu et al.ll2007l : IHughes et al.ll2008l . 120111 ) 

clouds {Xco ^ 10""^) is also applicable in the disk, the molecular density normalization parameter 
Mco implies a modest gas mass, 2.8t.i'l ^ IO^^Mq, validating a posteriori our assumption that 
self-gravity is negligible {M^/M^ ^ 10~^ <C 1). In any case, our simple treatment of the disk 
structure is really only used as a means to an end: the focus here is to place a firm constraint on 
the key parameter that determines the behavior of the disk velocity field - the central stellar mass, 
M*. 

We infer a total stellar mass of = 1.75^0 O6^0 V4046 Sgr spectroscopic binary, a 

formally precise constraint with a relative uncertainty of 3-5%. Figure S] demonstrates that this 
estimate of is essentially independent of the disk structure parameters in the model. The sole 
relevant degeneracy is with the disk inclination, which we infer to be id = 33? 5^54 relative 
precision of ^2-5%). Within t he quoted uncert a inties , these {M^, id} values are consistent with 
the previous determinations by [Rodriguez et al.l (I2OIOI ) , who used a simple g^id search and an 
altogether different underlying disk structure model. These {M*, id} measurements based on the 
circumstellar gas disk kinematics can be directly compared with the joint (degenerate) constraints 
on {M^, z*} imposed by RV measurem ents of the spectr o scopic binary itself . Although various 



group s have studied this SB2 system (lOuast et all [20001; IStempels fc Gahml l2004l : iDonati et al 



20111), by far the most robust constraints on the stellar parameters come from the l ong-term 
extens ive, and high resolution optical spectroscopic monitoring campaign conducted by IStempels 
(I2OI2I ). In that work, the RV data are found to be best explained with a binary that has a total 
mass M^{smuf = 0.2923 ± O.OOOTM©. Thejnset inFi gure [4] confirms that our disk kinematics 
constraints and the RV constraints by IStempelsl (l2012l ) coincide in the stellar mass— inclination 
plane: the two completely independent methods find the same results well within their formal 1 a 
uncertainties. Adopting our best estimate of and propagating the relevant uncertainties, the 
RV constraints s uggest that = 33?42 it 0.01. Therefore, we quantitatively confirm the finding of 
Rodriguez et al.l : the V4046 Sgr spectroscopic binary and its associated, large-scale circumbinary 
disk are co-planar within \id — i*| ^ 0.1 =t 1° across ^4 orders of magnitude in radial scale. 



- 12 - 



60 
54 
48 
42 

0.16 

0.14 

2.0 

1.8 

1.6 
36 

34 
32 

120 
110 

0.65 
0.60 



V 




# A 



(8? 




# 






35 



34 



33 



32 





this paper 
CO disk rotation 




Stempels 2012 
RV monitoring (SB2) 



1.6 



1.7 1.8 
[Mo] 





1.9 






4 8 42 48 54 60 0.14 0.16 1.6 1.8 2.032 34 36 110 120 0.60 0.65 



M,, [10-« Mo] 



Fig. 4. — Marginalized pararameters of tlie sampled posterior distrubution. The 68% confidence 
contour is given by the red contour with each binned distribution shown by the linear greyscale. 



5. Discussion 



We have presented spatially and spectrally resolved Submillimeter Array observations of the 
^^CO J=2— 1 line emission from the isolated, nearby, and gas-rich circumstellar disk that orbits 
the close, pre-main sequence double-lined spectroscopic binary V4046 Sgr. Adopting a simple 
parametric model for the structure of a flared Keplerian disk, we employ a Monte Carlo Markov 
Chain technique to infer the properties of the disk velocity field from these data, and thereby 
provide a robust statistical estimate for the total mass of the central binary. We find that these 
CO line data are best reproduced for a binary mass = I.TS^qq^ Mq and a disk viewing angle 
inclined by = 33? 5^^ 4 f roni face-on. Tho se values are in excellent agreement with the completely 
independent inferences of IStempelsl (120121 ). made from their extensive radial velocity monitoring 
campaign of the spectroscopic binary itself. The orbital planes of the binary (a = 0.045 AU) and its 
associated circumbinary disk (on radial scales out to ^400 AU) are aligned within ^0.1±1°. These 
results demonstrate that, despite its complexity, the disk-based kinematic method for estimating 
the masses of young stars is both precise (at the level of a few percent) and accurate in an absolute 
sense, as verified here by an entirely independent dynamical constraint. 



The disk-based dynamical estimate of can be combined with the RV constraints of lStempels 



(|2012l ) and other information from the literature to assess the predictions of pre-MS stellar evolution 
models. Coup ling our estimate with the stellar mass ratio (q) and mass function {M^(sini^)^) 
determined by Stempelsl . we can derive the masses of the individual stellar components in the V4046 
Sgr binary, 0.90 ±0.05 and 0.85 ± 0.04 Mq, and their orbital inclination, = 33?42±0.0li Assum- 
ing that the stellar rotation axes are aligned with the binary or bital axis (a spin-orbit alignment 
presu mably induced by tides in this close, circularized system; e.g., Zahnlll977 : Hutlll98ll:lMelo et al 



I2OO1I ). the radial and rotational velocities measured for each stellar co mponent by Stempe ^ 



can 



be used to determine the stellar radii, 1.25 ± 0.04 and 1.21 it 0.04 R©. iStempels fc Gahml (120041 ) 
inferred spectral types of K5 and K7 in this system, corresponding to effective temp e ratur es of 
4350 ± 240 and 4060 =b 210 K for the standard conversion advocated by ISchmidt-Kalerl (119821 ): an 
ambiguity of one spectral subclass has been assumed to estimate the temperature uncertainties. 
Those measurements imply that the individual stellar luminosities are 0.50 ±0.11 and 0.36 ±0.08 
(a total luminosity of 0.86 ± 0.14 L©). 

Figure [5] shows the composite optical/near-infrared spectrum expected from the V4046 Sgr 
binary, given the effective temperatures, stellar radii, and masses derived above from the com- 
bination of RV data and our disk-based measurements (black curve^ with uncertainty in shaded 
gray). This spectral energ y distribution (SED) prediction was generated from an interpolated 



grid of iLejeune et al.l (119971 ) model spectra, assuming the appropriate stellar gravities (log g 
4.20 ± 0.02 in both cases), negligible interstellar reddening {Ay — 0), and our adopted d — 
73 pc. Representative broadband photometric data are marked as the datapoints in Figure [5l con- 



^Here and throughout, we conservatively adopt the larger of all asymmetric uncertainties for simplicity. 
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Fig. 5. — The composite optical/near-infrared broadband spectrum observed for the V4046 Sgr 
binary (see text for references and discussion of uncertainties), along with a model spectrum pre- 
dicted from a combinati on of the CO ro tation curve analysis in §4 and the optical spectroscopic 
monitoring campaign of IStempelsl (l2012l ) (black curve, with uncertainties shaded in gray). The 
red curve corresonds to that prediction scaled down by ^30% in luminosity, to best reproduce 
the observations. The data and the predictions are marginally discrepant (1.4 a) in terms of their 
normalization; potential causes of that minor difference are discussed in the text. 



structed from the weight ed mean magnitud es of iHutchinson et al.l (ll990l ). IStrassmeier et al.l (119931) . 
Jensen fc MathieJ fll997l ). and the 2MASS flskriitskie et albood l and DENIS flEpchtein et al.lll997l l 
surveys: uncertainties were taken as the quadrature sum of the standard deviation of the weighted 
mean and the maximal deviation from the weighted mean, in an effort to more faithfully represent 
potential uncertainties due to variability. A quick examination of Figure [5] demonstrates that the 
spectral morphology of the data and model prediction are a good match, although the normaliza- 
tions appear discrepant. Scaling the (best-fit) predicted spectrum down by ^30±15% provides a 
good match to the observed spectrum; a composite luminosity of 0.60 ±0.13 Lq is more appropriate 
(red curve). Therefore, the predicted and observed composite binary luminosities are marginally 
(^1.4 a) different. 

Although formally the conflict in these luminosity values is not statistically significant, it is 
still interesting to consider some potential paths for a better reconciliation of all the data. Perhaps 
the most straightforward means of doing so is to modify the assumed distance to the system: taken 
at face value, a ^30% increase to d ^ 95 pc implies an inherently more luminous pair of stars 
with a spectrum that would be in good agreement with the observations. However, that same 
re-normalization impacts the total stellar mass inferred from the CO data, with a linear scaling 
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to ^ 2.3 M0 that introduces a substantial {^2 a) discrepancy in the {M*, z}-plane between 
the disk-kinematics and RV techniques for estimating stehar parameters. Although it hardly seems 
worthwhile to trade one marginal discrepancy for another (which is technically less marginal) , there 
is no a priori reason that the orbital planes of the disk and binary need to be aligned with the 
precision inferred in §4. Unfortunately, little guidance is pr ovided in the form of an uncertainty on 
the kinematic parallax measurement of Torres et al. fcood ). In any case, discrepanc ies in effective 



temperatures and luminosities are not necessarily uncommon f or active young; stars f|Stassun et al. 



20121), and particularly for those in close binary systems (e.g., iGomez Maqueo Chew et al.ll2012l ). 



Ultimately, dwelling on a marginal lumi nosity discrepancy i s not well justified, especially given 



the independent distance estimate from the iTorres et al.l (120061 ) study. In the following, we adopt 



the stellar parameters inferred from the joint constraints of the CO disk spatio-kinematics and the 
optical RV studies of Stempels. assuming = 73 pc, but use individual stellar luminosities scaled 
down by 30% to best match the observed optical and near-infrared spectrum: 0.35 z b 0.10 L^t^ and 



0.25 =b 0.08 L0 (note that each component is slightly less luminous than reported by IPonati et al. 



20111). With those results, we can make comparisons with the predictions of pre-MS evolution 



models to estimate the ages (t*) and masses of the individ ual stars in the V4046 S gr bin ary. To 
accomplish that, we follow the Bayesian methodology of *j0rgense n fc Lindegren fcood ) . which 



employs a finely-interpolated grid of pre-MS models in the HR diagram to estimate the probability 
distributions of stellar mass and age given the measured values (and uncertainties) of effective 
temperatures and luminosities (see also lCennaro et al.ll2012h . If we define the observables as {x, y} 
= {logT, logL} with associated uncertainties {ax, Cy}^ and the pre-MS model predictions {x(M*, 
t*), y(M^, ^*)}, we can write the likelihood function as a multivariate Gaussian, 

1 ^ (x-x)^ ^ {y-yf 
2 



C{x,y\x,y) 



1 



■ exp 



(12) 



The best-fit {M*, t*} are then directly determined from the {x, that maximize the likelihood, 
with uncertainties that can be calculated from the shape of the likelihood distribution. This 
procedure was conducted for four different pre-MS evolution models, assuming solar composition, 
a fractional deuterium abundance of 2 x 10~^, and a convecti v e mix ing; p arameter a ^ 1.7-2.0: 



D'Antona fc Mazzitellil (119981 ). iBaraffe et al.1 (119981 ). ISiess et all (12000), and iTo gnelli et al.1 (1201 ih . 



The results of these comparisons are presented in Table [2l Figure [6] shows the V4046 Sgr stellar 
properties that were inferred from each of the four reference sets of pre-MS evolutionary model 
tracks in the HR diagram. All of the evolutionary models make predictions for the V4046 Sgr 
stellar masses that are remarkably consistent with each other and the dynamical masses inferred in 
§4. The stellar masses that formally maximize the likelihood in Equation 12 are systematically ^2- 
10% below the best-fit dynamical masses, an un der-prediction typical of pre - MS models regardles s 
of the dynamical method used to estimate (^ Hillenbrand fc White 2004; Gennaro et al. 2012 ). 
However, those modest discrepancies are not statistically significant, given the uncertainties on {L, 
T} and the dynamical masses. An examination of Figure [6] demonstrates that the V4046 Sgr stars 
are nearly evolved off the Hayashi track, having presumably developed radiative cores as expected 
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Fig. 6. — (top) The V4046 Sgr stellar properties (pri mary in red, secondary i n blue) c ompared with 



the pre-MS evolution ary models of (fro m left to rigMllD'Antona fc Mazzitellil (DM98) , iBarafFe et al 



(BCAH98), ISiess et all (SDFOO), and iTognelli et all (1201 ih (TPDll) in the HR diagram. In each 
case, the mass tracks for = {0.7, 0.8, 0.9, 1.0} Mq are shown as gray solid curves, with isochrones 
for = {1, 10, 100} Myr denoted as gray dashed curves, (bottom) The 68% confidence intervals for 
the posterior likelihood distribution in the |M^ , t^j-plane derived froni the H R diagram for each 
model type, using the Bayesian methodology of |j0rgensen fc LindegrenI (120051 ) . The vertical lines 
indicate the individual stellar masses as determined by combining the tot al mass from th e disk 
analysis in §4 and the mass ratio derived from the spectroscopic analysis of IStempelsl (l2012l ). 



for their solar-like masses. The clear implication of their location in the HR diagram is that the 
V4046 S^T b i nary i s com paratively old i oi a p re-MS system, as has been previously reported by 



IDS 



Donati et al.l (120111 ) and iKastner et aLl (120111 ). The models used here suggest a large range of 
ages are plausible, from ^5-30 Myr, and confirm that the binary components are coeval. If we 
include a Gaussian prior representing the dynamical masses into the Bayesian analysis of the HR 
diagrams described above (labeled as '+dyn' in Table Ej), we can infer a smaller range of acceptable 
"dynamical" ages from each set of pre-MS models. The different pre-MS model predictions with 
these dynamica l priors are shown t ogether, and averaged (weighted by the posterior probability 
for each model; iHoeting et al.lll999l ). in Figure [2 the averaged results suggest ages of 12tY and 
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Fig. 7. — The 68% confidence intervals for tlie posterior likelihood distribution in the {M*, t^}- 
plane derived from the HR diagram for each model type (see key in right panel) overlaid together on 
the individual V4046 Sgr components {left primary; right secondary), now assuming the dynamical 
masses determined in §4 as Gaussian priors (the '+dyn' models references in Table Ej). The filled 
shaded region corresponds to the averaged likelihoods from the four different pre-MS models. The 
dynamical masses and their uncertainties are marked as a hatched, shaded vertical column. 



Myr for V4046 Sgr A and B, respectively. The corresponding coeval age is ISlg Myr, in good 



agreement with the age constraints from the p utative far-flung: com panion(s) V4046 Sgr C[D] (at 
a separation of ^12,350 AU, or 2f8 on the sky; 



Kastner et al. 



20111). 



Torres et al.l (120061 ) suggested that the observed space motion of V4046 Sgr is consistent with 
a high probability of membership in the /3 Pic moving group, a widespre ad and local (d ^ 34 zb 



21 pc) stellar association with an age range e stimated to be ^8-20 Myr (e.g. JZuckerman et al.ll2001 



Zuckerman fc Song l2004l : iTorres et al.ll2008l ). The kinematic parallax calculated by Torres et al. 
{d ^ 73 pc) and component ages estimated here for the V4046 Sgr binary are certainly supportive of 
that conclusion. If that is the case, it is worth noting the uniqueness of the V4046 Sgr circumstellar 
environment: there are no other /3 Pic moving group members that are known to host such a 
massive, long-lived, gas-rich accretion disk. Regardless of its original association, the proximity 
and advanced age of V4046 Sgr make for a remarkable case study in both the long-term evolution 
of protoplanetary disk structure and the fundamental properties of pre-MS binary stars. 

Although V4046 Sgr is a particularly striking example of the methodology behind disk-based 
dynamical estimates of stellar masses, we anticipate that these techniques will find considerably 
more use in the near future as the Atacama Large Millimeter Array (ALMA) project is completed. 
Even with vastly improved data quality, this simplified model to extract the gas disk rotation curve 
from such interferometric observations (see §3) remains complex. However, we have demonstrated 
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that the method is robust, accurate, and precise, by using independent dynamical constraints 
from the V4046 Sgr spectroscopic binary radial velocity m onitoring results. In practice, these 
constraints and others like them (notably for UZ Tau E; see lSimon et alj|2000l : IPrato et al.ll2002l ) 
effectively serve as a check on the absolute calibration of the modeling technique described in §3. 
The consistency between the independent dynamical constraints on the V4046 Sgr stellar masses 
validates the application of this procedure to dynamically "weigh" single, isolated pre-MS stars 
based on their gas disk kinematics. Along with an investment in accurate stellar luminosity and 
temperature measurements, ALMA observations of the molecular gas kinematics in circumstellar 
disks will usher in a new era of precision in the fundamental astrophysical properties of young stars. 



6. Summary 

We have presented sensitive, high-resolution SMA observations of the ^^CO J=2— 1 line emis- 
sion from the massive, gas-rich disk orbiting the double-lined spectroscopic binary star V4046 Sgr. 
Using simple radiation transfer calculations for a disk structure model with a Keplerian velocity 
field, we fit the observed spectral visibilities using a stochastic model optimization technique that 
simultaneously infers model parameter values and their uncertanties. Our specific focus has been 
on the key model parameters that describe the disk velocity field, with a goal of placing a firm, 
dynamical constraint on the mass of the central binary. The key conclusions of our analysis are: 

1. From modeling the CO line emitted by the circumbinary disk, we infer that the total stellar 
mass of the V4046 Sg r binary is = assuming the kinematic parallax distance 



of 73 pc estimated by iTorres et al.l (120061 ). That measurement is in excellent agreement with 



the independent dynamical constraints imposed from the analysis of an opti cal spectroscopi c 



monitoring campaign of the radial velocity variations from the binary itself (IStempelsll2012h . 



The mutual consistency of these distinct dynamical constraints on the stellar masses verifies 
that the mass determination based on the velocity field of the disk gas is accurate in an 
absolute sense, as well as remarkably precise (with a 3-5% formal uncertainty on M*). 

That same combination of constraints from millimeter and optical spectroscopic measure- 
ments confirm that the orbital planes of the stars and their accompanying circumbinary disk 
are co-planar, with inferred inclination angles {id = 33?5l5'4 ^* — 33?4±0.01) that differ 
by 0.1 ± 1° over roughly four decades in radius, from ^0.04 to 400 AU. 

The inferred component masses of the V4046 Sgr binary are in uniformly good agreement 
with a variety of pre-MS evolution model predictions. When combined with our dynamical 
constraints, those same models confirm the coevality of the binary components, and suggest 
an average age for the system of 13^3 Myr. Therefore, V4046 Sgr hosts one of the oldest and 
nearest gas-rich primordial accretion disks currently known. 
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Table 1. Disk Model Parameters 



parameter units estimate 





[10-6 Mq] 




rc 


[AU] 


45^1 


q 




0.63 ±0.01 




[K] 


115 ±5 


e 


[km s"-*-] 


0.14 ±0.01 


id 


[°] 


33.5l?-I 




[Mq] 


-■■•'^-.oe 



Note. — See §3 for a description of 
the model parameters. The uncertain- 
ties for each parameter are determined 
from the width of the posterior distribu- 
tion that encapsulates 68% of the mod- 
els around the peak (equivalent to la 
for Gaussian random variables. 



Table 2. Comparison with pre-MS Evolution Models 



V4046 Sgr A V4046 Sgr B 

[Mq] U [Myr] ~~ U 

Models HR +dyn HR +dyn HR +dyn HR +dyn 



DM98 
BCAH98 
SDFOO 
TPDll 



+0.06 
-0.10 
+0.09 
0.12 
+0.13 
0.06 

0.86tH? 



0.84 
0.89 
0.83 



0.88 
0.90 



+0.03 
-0.05 
+0.04 



-0.06 

0.881°:°^ 
0.891°:°^ 



6if 

lOlf 
12lf 

9+35 



9: 
18: 



-10 
-2 
-10 



-9 

14+10 



0.76 
0.80 



+0.06 
-0.09 
+0.07 
-0.09 

0.75tH^ 
0.78tHI 



0.82 
0.84 



+0.03 
-0.05 
+0.04 
-0.05 

n Oq+0.04 

'-'•^'^-0.05 
'-'•^^-0.04 



8lf 
9tf 



(mean) 0.84 



+0.11 
-0.08 



0.89 



+0.04 
-0.05 



1+43 



12 



-17 



0.77 



+0.08 
-0.11 



0.83 



+0.04 
-0.05 



Note. — Stellar mass and age estimates from the four different pre-MS evolution models 
are calculated using the HR observables {logL, logT}, following the methods outlined in 
§5. We report results for both a uniform (HR) and Gaussian (+dyn) prior on the mass, 
with the latter based on the dynamical constra ints. The mean value s are computed for all 
the models, following the technique outhned bv iHoeting et al.l (|l999l ). 
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